1 Introduction

With the development of advanced 3D CAD tools in the last decades, architects are now able to incorporate complex freeform surfaces into their designs. This new capabilities have given birth to new architectural possibilities and, more importantly, new needs and challenges for the industry that can often be solved by the means of Differential Geometry (Carmo 2016) and, more specifically, Discrete Differential Geometry (Keenan Crane 2013), which studies the discrete counterparts of Differential Geometry. This new fusion between architecture and mathematics gave birth to the field known as Architectural Geometry (Pottmann, Schiftner, and Wallner 2008; Pottmann 2010), which has been further developed during the past years.

Given the scale architectural projects, this smooth freeform surfaces cannot be built as designed in 3D CAD applications and, therefore, must be divided in smaller, simpler parts that could be easily manufacturable. This process is called Rationalization, and is not a trivial task; since the Rationalization of a surface implies a loss of smoothness in the resulting shape, it must be done in such a way that minimizes this loss. The selected shape of the generated panels also plays a key role when approaching this problems, since it is clear the final appearance of the design will vary greatly depending on wether the panels are triangular, quad, hexagonal or long strips. Constructibility of the selected design (panel planarity, curvature, error margins) is of great concern in rationalization methods for architectural purposes, looking for robust and precise algorithms, as opposed to related Computer Graphics methods, whose main focus is usually related with execution times & speed.

In this publication we center on the generation of Geodesic Panels on Freeform Surfaces, meaning, seamlessly covering any given surface with a pattern of planks; thin, long, straight panels, which must be developable and rectangular (or quasi-rectangular) in shape and must approximate the surface by being bent only on their weak axis. We will introduce the basic geometric properties behind this concept and the algorithmic implementation of several of the methods. We will also explain some concepts developed for optimizing any given shape towards developability.

2 Background

Examples of the application of Architectural Geometry are abundant in today’s architecture. Architect Frank Ghery is known for his work ground-braking work on the design of architectural shapes based purely on developable patches (Stein, Grinspun, and Crane 2018), which can be easily covered by a variety of panel patterns. Many of his buildings have become iconic landmarks of the cities they are built in

Covering surfaces with long rectangular, or quasi-rectangular panels of equal width, is a topic widely explored in traditional boat building, and has also been recently applied in several architectural projects, such as Toyo Ito’s Minna No Mori (Japan), or the interior cladding at Frank Ghery’s Burj Khalifa (Dubai) (Meredith and Kotronis 2013).

NOX Architects, an award winning architecture firm, is famous for building complex shapes using planks. Their design methodology included the creation of very precise models, and manually covering these models with paper strips until a final cladding solution could be achieved fig. 5

Some recently finished buildings would have also benefitted from the optimization of paneling and/or beam layout like the Yeoju Golf project by Shigeru Ban, which involved the construction of complex laminated wood shapes fig. 6

Geodesic curves on manifolds and their calculation is also important for this topic (V. Surazhsky et al. 2005; Cheng et al. 2016,Polthier and Schmies (1998)) as well as the concept of developability of surfaces and their discrete counterparts (see Carmo 2016, and more recently Stein, Grinspun, and Crane (2018)). In this publication, all of the geodesics will be generated with the initial value problem.

The concept of mesh developability has also been widely studied in the past, although strictly speaking, it will not be the object of this publication. Previous attempts at mesh developability can be viewed in fig. 7. The latest and more robust implementation was developed by Stein, Grinspun, and Crane (2018) and is based on reducing the angle defect on each vertex on a mesh, by forcing the normals of the faces adjacent to each vertex to be parallel; this causes that each vertex on the given mesh becomes either flat or a hinge fig. 8 (a), effectively concentrating the gaussian curvature of the original mesh on the hinges. This procedure provides very smooth developable results on any given mesh, but most of the time, the resulting developable will exhibit creases and sharp edges in some areas, which are heavily influenced by the topology of the starting mesh. Also, the resulting flat configuration of the mesh usually contains self-intersections and complex boundary shapes that would increase complexity on the construction and assembly stages in an architectural scale project.

3 Geometric properties

3.1 Geodesic curves

In differential geometry, a geodesic curve is the generalization of a straight line into curved spaces. It is usually understood as the shortest path between two points on a given surface, although this may not always be the case.

Also, in the presence of an affine connection, a geodesic is defined to be a curve whose tangent vectors remain parallel if they are transported along it. We will explore the notion of vector parallel transport in the following sections.

For triangle meshes, shortest polylines cross edges at equal angles.

Finding the truly shortest geodesic paths requires the computation of distance fields (see Carmo 2016; Kimmel and Sethian 1998)

Algorithmic ways of generating geodesics

The computation of geodesics on smooth surfaces is a classical topic, and can be reduced to two different solutions, depending on the initial conditions of the problem (Deng 2011):

Initial value problem

given a point \(p \in S\) and a vector \(v \in T_pS\), find a geodesic which is incident with \(p\), such that the tangent vector fo \(g\) at \(p\) is \(v\).

Boundary value problem

given two points \(p_1, p_2 \in S\) find a geodesic \(g\) connecting \(p_1\) and \(p_2\).

Both problems have different ways of being solved either numerically, graphically or by the means of simulations. The initial value problem can be solved using the concept of straightest geodesics (Polthier and Schmies 1998), whereas the boundary value problem has a very close relation with the computation of shortest paths between two points on a surface.

3.2 Developable surfaces

It is also important to introduce the concept of the developable surface, a special kind of surface that have substantial and variable normal curvature while guaranteeing zero gaussian curvature, and as such, this surfaces can be unrolled into a flat plane with no deformation of the surface. This surfaces have been an extremely important design element in the practice of known architect Frank Ghery and explained in (Shelden 2002).

There are several ways of generating a developable surface using a curve in space. The most simple way is the following.

1. Select the starting point of the curve
2. Obtain the perpendicular frame of the curve at the specified parameter.
3. Deconstruct the frame into it’s \(X,Y\) & \(Z\) vectors.
4. Select one of the vectors at all the sample points
5. Draw a line using the selected vector at the start point of the curve.
6. Create a surface using the rulings obtained in step 5.

Not all developable surfaces are of interest in an architectural sense. We must make a distiction between noisy and smooth developable surfaces, as can be seen in Fig. 10. We are spetially interested in generating smooth develpable surfaces, which can be easily manufactured, as opposed to the complex bending operations that manufacturing a noisy developable surface will require

3.3 Properties to aim for

We will first introduce the properties we aim for when generating geodesic panels:

Geodesic property

• Long Thin panels that bend about their weak axis
• Zero geodesic curvature
• Represent the shortest path between two points on a surface

Constant width property

• Panels whose original, unfolded shape is a rectangle.
• The only way this can happen is if the entire surface is developable.
• For all other surfaces:
  • Assuming no gaps between panels
  • Panels will not be exactly rectangular when unfolded
  • Requirement: Geodesic curves that guide the panels must have approximately constant distance from their neighborhood curves.

Developable (or ‘pure-bending’) property

• Bending panels on surfaces changes the distances in points only by a small amount so,
• A certain amount of twisting is also present in this applications.
  Some methods in this chapter do not take into account this property.

Objectives

Look for a system of geodesic curves that covers a freeform surface in a way that:

1. They have approximate constant distance with it’s neighbours.
2. This curves will serve as guiding curves for the panels.
3. The panels are to cover the surface with no overlap and only small gaps

4 Algorithmic strategies

In this section, we will introduce the different existent methods for generating geodesic curve patterns that will be used to model the final panels, using methods indicated in sec. 4.6.

4.1 Parallel Transport Method

This method, described in (Pottmann et al. 2010), allows for the generation of a system of geodesic curves where either the maximum distance or the minimum distance between adjacent points occurs at a prescribed location.

In differential geometry, the concept of parallel transport (see fig. 11) of a vector V along a curve S contained in a surface means moving that vector along S such that:

1. It remains tangent to the surface
2. It changes as little as possible in direction
3. It’s length remains unchanged

Procedure

Following this procedure, extremal distances between adjacent geodesics occur near the chosen curve. Meaning:

1. For surfaces of positive curvature, the parallel transport method will yield a 1-geodesic pattern on which the maximum distance between curves will be \(W\).
2. On the other hand, for surfaces of negative curvature, the method will yield a 1-geodesic pattern with \(W\) being the closest (or minimum) distance between them.

The placement of the first geodesic curve and the selection of the initial vector are not trivial tasks. For surfaces with high variations of surfaces, the results might be unpredictable and, as such, this method is only suitable for surfaces with nearly constant curvature. Other solutions might involve cutting the surface into patches of nearly-constant curvature, and applying the parallel transport method independently on each patch.

4.2 Evolution Method

Two main concepts are covered in this section, both proposed by (Pottmann et al. 2010): the first, what is called the evolution method, and a second method based on piecewise-geodesic curves, curves that are not geodesic, but are composed of several connected geodesic segments.

The evolution method

Starting from a source geodesic somewhere in the surface:

• Evolve a pattern of geodesics iteratively computing ‘next’ geodesics.

• ‘Next’ geodesics must fulfill the condition of being at approximately constant distance from its predecessor.

• ‘Next geodesics’ are as described in fig. 15

Distances between geodesics

1. No straight forward solution.
   1. Only for rotational surfaces (surfaces with evenly distributed meridian curves).
2. But a first-order approximation of this distance can be approximated:

The method proposed in Pottmann et al. (2010) is based on the calculation of Jacobi Fields, vector fields perpendicular to a geodesic that approximate the distance between one geodesic and the next in an infinitesimal distance.

Starting at time \(t=0\) with a geodesic curve \(g(s)\), parametrized by arc-length \(s\), and let it move within the surface.
A snapshot at time \(t=\varepsilon\) yields a geodesic \(g^+\) near \(g\).

\[g^+(s)=g(s)+\varepsilon\mathbf{v}(s)+\varepsilon^2(\ldots)\qquad(1)\]

The derivative vector field \(\mathbf v\) is called a Jacobi field. We may assume it is orthogonal to \(g(s)\) and it is expressed in terms of the geodesic tangent vector \(g'\) as:

\[\mathbf{v}(s)=w(s)\cdot{R_{\pi/2}(g'(s))},\quad\text{where}\;w''+Kw=0\qquad(2)\]

Since the distance between infinitesimally close geodesics are governed by Eq. 2, that equation also governs the width of a strip bounded by two geodesics at a small finite distance.

We propose a simpler approach to this problem, since the computation of the jacobi field can give very high margins of error and ‘next’ geodesic \(g_{i+1}\) will still need to be computed using the initial value problem. The modification is the following: we can substitute the length of the jacobi field for the user’s specified width \(W\) and find the best fit geodesic that fits this signed distances.

Obtaining the ‘next’ geodesic

Any given surface can be completely covered in this manner by recursively computing the next geodesic using the previous, given some margin of error.

Piecewise evolution method

We can modify the implementation of the evolution method to be able to control the distance between geodesics even in areas with a high rate of change in curvature. In order to prevent this unwanted width variations, instead of looking for the best overall geodesic, we will look for the best geodesic that fits the largest interval of \(h_i\)’s possible fig. 16, given a specified threshold \(\epsilon\). Once this interval is obtained, the current geodesic is split and a new geodesic will be generated. This process will be repeated until the curve crosses all \(h_i\)’s (some exceptions apply).

4.3 Level-Set Method

We can also compute geodesic curve patterns on a surface as a series of levels (or level-sets) of a scalar-valued function assigned on each vertex of a mesh. This concept was first introduced in Wallner et al. (2010) and is based on concepts developed in the geodesic active contour method.

The computation of level-sets on a triangle mesh is a fairly simple task. Since the values of the function are assigned per vertex, the line of a constant value of the function is obtained by linear interpolation on each of the mesh faces fig. 18. The gradient of the scalar-valued function also needs to be calculated, which results in a vector field which is constant on each face. Once the gradient has been computed, the geodesic curvature of a level set can be obtained as the divergence of the gradient vector field over each vertex:

\[K_g = div(\frac{\nabla \phi}{||\nabla \phi ||})\]

In order that the obtained levels approximate the surface, we can approach the problem as the minimization of a combination of error measurements:

\[F_{min} = F_k + \lambda F_\nabla + \nu F_w\]

where \(F_k\) measures the deviation of the scalar function geodesic property, \(F_\nabla\) is computed for regularization purposes and \(F_w\) measures the constant-width deviation; \(\lambda\) and \(\nu\) are user defined weighting values. \(F_k\) and \(F_w\) will only vanish at the same time in the case of developable surfaces. This values can be computed on a mesh (V,E,F) using the following formulas, explained with more detail in Pottmann et al. (2010):

\[F_k = \sum_{\mathbf{v}\in V}{\mathcal{A}(\mathbf{v})\bigg(div(\frac{\nabla \phi}{||\nabla \phi ||}(\mathbf{v}))\bigg)^2}\] \[F_\Delta = area(S)\sum_{\mathbf{v}\in V}{\mathcal{A}(\mathbf{v})\Delta \phi(\mathbf{v})^2}\] \[F_w = \sum_{f\in F}area(f)\bigg(||\nabla\phi(f)||- \frac{h}{w}\bigg)^2\]

The algorithm would be started with a set of randomly placed values on each vertex of the mesh fig. 19. The implementation used in this paper differs from the originally proposed in Pottmann et al. (2010), in which the Gaussian-Newton method is used for minimization and all high order derivatives are computed analytically. In our implementation, we explored the application of gradient-less minimization algorithms, in particular BOBYQA (Bound Optimization BY Quadratic Approximation) although this decision causes computation times to increase severely.

4.4 Geodesic Webs

We can also expand the Level-Set method explained in the previous section from one family of curves to several families of interconnected quasi-geodesic curves by simply computing several sets of scalar functions on each vertex of the surface.

This can be easily done by adding another equation to our minimization problem.

\[\nu \dot F_{angle} = \nu \dot \sum_{f\in F} \frac{area(f)}{area(S)}\langle \frac{\nabla\phi_i}{||\nabla\phi_i||}.R_{\frac{\pi}{2}-\alpha}(\frac{\nabla\phi_j}{||\nabla\phi_j||})\rangle^2\qquad(3)\]

where \(\nu\) is the weight given to \(F_{angle}\), \(R_{\frac{\pi}{2}-\alpha}\) means the rotation by the angle \(\frac{\pi}{2}-\alpha\) in the respective face.

4.5 Geodesic Vector fields

The level-set method described in Sec. 4.3 is not suitable for arbitrary surfaces, and therefore it must be adapted to achieve the desired result. In this section, we introduce the concept of Geodesic Vector-fields (Pottmann et al. 2010) to divide the surface into patches that could be easily covered by equal width panels using the level-set method.

A geodesic vector field on a surface \(S\) is geodesic if it consists of tangent vectors of a 1-parameter family of geodesics covering \(S\) (Pottmann et al. 2010).

Considering the surface \(S\) as a triangle mesh (V,E,F), the vector field is represented by unit vectors \(\mathbf v_f\) attached to the incenters of facef \(f \in F\).

\[\mathbf v_{f_2} = \mathbf v_{f_1} + J_{f_1} \dot (m_{f_2} - m_{f_1}) + r_e\qquad(4)\]

\[ where J_f = XXXX\]

The vector field \(\mathbf v = (\mathbf v_f)_{f \in F}\) is a geodesic vector field if it’s length is constant for all faces \(f\), and we can find a collection of coefficients \(\mathbf g = (\mathbf g_f)_{f \in F}\) with \(\mathbf g_f = (g_{1,f},g_{2,f},g_{3,f})\), such that eq. 4 holds with \(\mathbf r\)’s bellow some theshold.

The following equation quantifies how well \(\mathbf v\) satisfies the geodesic property:

\[ Q(\mathbf v) = min_{g\in R^{ 3|F| }}\qquad(5)\]

4.6 Panel modeling

In this section, we will discuss several ways to generate panels from curves lying on a given surface:

The Bi-Normal Method

The second method for defining panels, once an appropriate system of geodesics has been found on \(\Phi\), works directly with the geodesic curves.

Assume that a point \(P(t)\) traverses a geodesic \(s\) with unit speed, where \(t\) is the time parameter. For each time \(t\) there is:

• a velocity vector \(T(t)\)
• the normal vector \(N(t)\)
• a third vector \(B(t)\), the binormal vector

This makes \(T.N.B\) a moving orthogonal right-handed frame

Although this method generates mathematically correct developable surfaces, Pottmann et al. (2010) describes it is unclear that the boundaries of the panels generated with this technique will correspond to real constructive requirements for the design, so we concentrate our efforts on method explained in the following section.

Tangent-developable method

The notion of Conjugate tangents on smooth surfaces must be defined:

• Strictly related to the Dupin Indicatrix
• In negatively curved areas, the Dupin Indicatrix is an hyperbola whose asymptotic directions (A1, A2)
• Any parallelogram tangentially circumscribed to the indicatrix defines two conjugate tangents T and U.
• The asymptotic directions of the Dupin indicatrix are the diagonals of any such parallelogram.

Initial algorithm is as follows:

For all geodesics \(s_i\) in a given pattern:

1. Compute the tangent developable surfaces \(\rightarrow\Psi_i\)
2. Trim \(\Psi_i\) along the intersection curves with their respective neighbours.
3. Unfold the trimmed \(\Psi_i\), obtaining the panels in flat state.

Unfortunately, this method needs to be refined in order to work in practice because:

1. The rulings of tangent developables may behave in weird ways
2. The intersection of the neighboring \(\Psi_i\)’s is often ill-defined.

Therefore, the procedure was modified in the following way:

This adjustments to the algorithm allow for the modeling of panels that meet the requirements of developability and approximately constant width, although it must be noted computation times increase, as double the amount of geodesic curves need to be generated, and subsequently, the desired width function needs to be half the desired width of the panels. Fig. 24 shows the result of computing such panels and the subsequent gaps between the generated panels; these gaps need to be kept within a certain width in order to produce a successful watertight panelization of the original surface. Furthermore, material restrictions such as bending or torsion where not taken into account during the construction of the panels.

5 Quality Assessment

5.1 Panel gaps

Since we are generating the panels using the tangent developable method, it is inevitable that some gaps appear between panels, especially when there is a high degree of curvature in the area. We can take this into account when calculating the trim curves, and obtain the maximum gap values of the model. The gap width is strictly related to the thickness of the panel, so if the maximum width of the gaps exceeds a desired a specified distance, computing narrower panels under the same conditions should guarantee the watertightness of the solution.

5.2 Stress and strain in panels

The following section investigates the behavior of a rectangular strip of elastic material when it is bent to the shape of a ruled surface \(\Psi\) un such way that the central line \(m\) of the strip follows the ‘middle geodesic’ \(g\) in \(\Psi\). This applies to both methods defining panels.fig. 25

• Lines parallel to \(m\) at distance \(d/2\) are not only bent but also stretched.
• If you introduce the radius of Gaussian curvature \(\rho = 1 / \sqrt{|K|}\)…
• …the relative increment in length of the strip is:

\[\varepsilon=\frac{1}{2}(d/2\rho)^2+\cdots\qquad(6)\]

Tensile stress

Tensile stress is estimated by expressing stress as \(\sigma = E\varepsilon\), which yields

\[ d/2\rho\leq C,\quad{with}\quad{C=\sqrt{\sigma _{max}/E}}\qquad(7)\]

where \(\rho_{max}\) is the maximum admissible stress and is Young’s modulus. Since this calculation is an approximation, a safety factor must be used when choosing \(\rho_{max}\).

\(C\) is a material constant, from which we can obtain an upper width bound for the panels by \(d_{max} = 2\rho_{min}C\), which will be the maximum admissible width for that design configuration.

Bending and shear stress

Bending (\(\sigma\)) and shear stresses (\(\tau\)) for panels with a thin rectangular cross-section depend on the panel thickness \(h\), but not on the panel width \(d\) if \(h/d \ll 1\), and maximum values occur on the outer surface of the panel (Wallner et al. 2010).

These values depend on the curvature \(\kappa\) of the central geodesic and the rate of torsion \(\theta\) of the panel. We have:

\[\sigma = E \kappa h/2 \quad and\quad \tau = hG\theta\qquad(8)\]

It is important to note that \(\tau\), measured by arc per meter, does not exceed \(\sqrt{|K|} = 1/\rho\), where the maximum value occurs in case the central geodesic’s tangent happens to be an asymptotic direction of the panel surface (Carmo 2016).

It is standard procedure to combine all stresses (tension, shear, bending) and use this information for checking panel admissibility.

Also, as one might conclude straight away, panels generated using the tangent developable method will experience less shear than the ones created by the ‘binormal’ method.

From here we can also check the admissibility of the computed models using the formulas described above.

Since the examples in this publication where created for illustration purposes, many of the models will not be admissible for fabrication. However, once a max panel width is calculated, if the design is not admissible, one can calculate narrower panels using the maximum calculated width in order to have an admissible model.

6 Conclusions

TEMPORARY

1. Parallel transport method and evolution method is only useful for simple surfaces of constant curvature, where the behavior of the generated geodesics can be easily predicted.
2. Piecewise evolution method is useful for any kind of surface, but there is a lack of control over the breakpoint positions, which might be undesirable depending on the design intention of the architect.
3. Level-set method, coupled with the geodesic vector fields, is a perfect solution to enable user interaction during the pattern design process, although in our implementation, we did not achieve sufficiently low computation times for this to be a possibility. It will also be necessary to check the deviation of the resulting curves’ geodesic property, since depending on the reference surface and the user specified weights, the generated pattern can deviate from true geodesics.

References